;******************************************
; plot regions with climate regime change 
;********************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;********************************************
 begin

 dirPath = ".../NFood_replication/data/" ; change to local path of your directory

 nmon = 12;
 syr = 2006;
 fyr = 2099;
 csyr = "" + syr
 cfyr = "" + fyr
 nyr = fyr - syr + 1
 dmis = 9.96921e+36

 ;case = (/"RCP45","RCP85","RCP85_SAI","RCP85_MSB","RCP85_CCT"/);
 ;case = (/"RCP85", "RCP45", "RCP85_SAI","RCP85_MSB","RCP85_CCT"/);
 case = (/"RCP85", "RCP85_SAI","RCP85_MSB","RCP85_CCT", "RCP45_85LU"/);
 ncase = dimsizes(case);

 ;8 major crops in CLM5: (NOT ACTIVE 15=C3 Generic crop),17=Temperate Corn, 19=Spring Wheat, (NOT ACTIVE 21=Winter Wheat),
 ;23=Temperate Soybean, 41=Cotton, 61=Rice, 67=Sugarcane, 75=Tropical Corn, 77=Tropical Soybean)
 ;cft_c3_crop = 1; cft_temperate_corn = 3; cft_spring_wheat = 5 ... cft sequence in crop landunit
 cftname = (/"C3 Generic crop", "Temperate Corn", "Spring Wheat", "Winter Wheat", "Temperate Soybean", \
             "Cotton", "Rice", "Sugarcane", "Tropical Corn", "Tropical Soybean" /)
 pfts = (/ 17, 19, 23, 41, 61, 67, 75, 77/)      ; crop index among all pfts
 cfts = (/ 17, 19, 23, 41, 61, 67, 75, 77/)-15   ; cft data index
 ncft = dimsizes(cfts)

;read dimension
 fdim = addfile(dirPath+"/land_cover/RCP85.1deg.landarea.nc", "r");
 lat = fdim->lat
 lon = fdim->lon
 nlat = dimsizes(lat);
 nlon = dimsizes(lon);
 dims = (/nlat,nlon/);
;print(lon)

;read data 
 ydbeg = new((/ncase,nlat,nlon/),"float"); Control period 2006-2015 for RCP85, RCP45
 ;dd2010@_FillValue = default_fillvalue("float")
 ydmid = new((/ncase,nlat,nlon/),"float");
 yd2075s = new((/ncase,nlat,nlon/),"float");
 ;yc2050 = new((/ncase,nlat,nlon/),"float");
 ;yc2095 = new((/ncase,ncft,nlat,nlon/),"float");
 ycmid = new((/ncase,nlat,nlon/),"float");
 yc2075s = new((/ncase,nlat,nlon/),"float");
 
 do kk=0,ncase-1
   fi := dirPath+"/yield_data/"+case(kk)+"_yield_avg_tC-ha_nyr.ncft.lat.lon" ;  [nyr,ncft,nlat,nlon]
   dd := fbindirread(fi, 0, (/nyr,ncft,nlat,nlon/), "float" ) ;read binary; nyr=94 (2006-2099)
   dd@_FillValue = default_fillvalue("float")
   ;Or wtave six crops first and then diff
   f2 := dirPath+"/yield_data/"+case(kk)+"_crop_area_ha_nyr.ncft.lat.lon"
   area := fbindirread(f2, 0, (/nyr,ncft,nlat,nlon/), "float" ) ;
   area@_FillValue = default_fillvalue("float") ; important to define NA, otherwise tp/ta will have devision by 0 error
   tp = dim_sum_n_Wrap(dd*area, 1); tot production per grid, tC
   ta = dim_sum_n_Wrap(area,1); tot crop area per grid, ha

   ;to avoid division by zero, use where() to mask out small areas for all scenarios
   yd = tp/where(ta.ne.0, ta, ta@_FillValue) ; tC/ha, wt ave yield (/nyr,nlat,nlon/)
   yd = where(ta.gt.10000, yd, ta@_FillValue) ;mask small crop cells <10000ha or 1% of grid cell
 
   ;every 20-year slice
   yd0 := yd(0:19,:,:)   ;CNTL time period: 2006-2025
   yd1 := yd(34:53,:,:) ;future time slice 1: 2040-2059
   yd2 := yd(69:93,:,:) ;future time slice 2: 2075-2099 
   ydbeg(kk,:,:) = dim_avg_n_Wrap(yd0,0) ;
   ydmid(kk,:,:) = dim_avg_n_Wrap(yd1,0) ;20-yr avg along the first dim
   yd2075s(kk,:,:) = dim_avg_n_Wrap(yd2,0) ;
   ycmid(kk,:,:) = ydmid(kk,:,:)-ydmid(0,:,:) ; yeild change relative to RCP85
   yc2075s(kk,:,:) = yd2075s(kk,:,:)-yd2075s(0,:,:)
   
 end do

 totyd = yd2075s/0.45*0.85 ;tC/ha to tonnes grain/ha yield for each case
 totyc = yc2075s/0.45*0.85 ;yield change relative to RCP85 or RCP45_85LU
 ;yc_rcp85 = (yd2095(0,:,:)-yd2050(0,:,:))/0.45*0.85 ; yield loss of rcp85 in 2090-2099 relative to 2046-2055

 ;printVarSummary(totyc) 
 print("range of yield change per grid: ")
 printMinMax(totyc,0) ;  
 ;min=-2.30877   max=2.06254
 ;print("mean="+avg(totyc))
 ;print("range of yield per grid: ")
 ;printMinMax(totyd,0)
 ;print("yield loss of rcp85 near end of 21cent: ")
 ;printMinMax(yc_rcp85,0)

 ;show more stats
 opt = True
 opt@PrintStat = True
 ycstat = stat_dispersion(totyc, opt) ;detailed statistics
;(0)      [23]     Lower 1.0%=-5.6576
;(0)      [24]     Lower 5.0%=-0.71055
;(0)      [25]     Upper 5.0%=1.92275
;(0)      [26]     Upper 1.0%=5.09478
 ;set color range between lower/upper 5% and min/max
 ;1deg*1deg = 1Mha; annual yield 1~10 tonnes/ha biomss
 ;set color range to -2 ~ +2 tonnes/ha yield

 totyc!0="case"
 totyc!1="lat"
 totyc&lat = lat
 totyc!2="lon"
 totyc&lon = lon
 totyc1 = lonFlip(totyc); 0~360 to -180~180

 copy_VarCoords(totyc, totyd);
 totyd1 = lonFlip(totyd); 0~360 to -180~180 flip together with data
 ;print(totyc1&lon) 
 ;copy_VarCoords(totyc(0,:,:), yc_rcp85);
 ;yc_rcp85 := lonFlip(yc_rcp85); 0~360 to -180~180 flip together with data


;plot
 imagePath = dirPath+"/../figures/supplementary"
 imagename = imagePath+"/Fig.5"
 wks_type = "pdf"        ; or "svg"
; wks_type = "png"
; wks_type@wkWidth = 1600 ;change resolution from default 1024x1024
; wks_type@wkHeight = 1600
 ;wks_type@wkOrientation = "landscape"; only works with ps or pdf
 wks_type@wkPaperHeightF = 8.5 ; default size 8.5x11 
 wks_type@wkPaperWidthF  = 11
 wks_type@wkDeviceLowerX = 0 
 wks_type@wkDeviceLowerY = 0
 wks_type@wkDeviceUpperX = 792
 wks_type@wkDeviceUpperY = 612
 wks1 = gsn_open_wks(wks_type,imagename)
 ;wks1 = gsn_open_wks("pdf",imagename)


 plot1 = new(ncase-1,graphic)

 res = True
 res@gsnDraw              = False
 res@gsnFrame             = False
 res@lbLabelBarOn         = False                 ;turn off individual label bars 

 res@cnLineLabelsOn       = False                 ; turn off line labels
 res@cnFillOn             = True                  ; turn on color fill
 ;res@cnFillPalette        = cmap                 ; cnFillPalette spans the color map by default
 ;res@gsnSpreadColors      = True                 ; this based on an old color scheme that was deprecated
 res@lbLabelAutoStride    = True                  ; optimal labels
 res@cnLinesOn            = False
 res@cnFillMode = "RasterFill"
 ;res@cnFillMode = "CellFill"

 res@gsnStringFontHeightF  = 0.025
 res@tmXBLabelFontHeightF  = 0.020
 res@tmYLLabelFontHeightF  = 0.020
 res@tmBorderThicknessF = 1
 res@tmXBMajorThicknessF = 1
 res@tmYLMajorThicknessF = 1

 res@mpCenterLonF = 0.
 res@mpMinLatF = -60.
 ;res@mpMaxLatF =  60.

 res@mpShapeMode  = "FreeAspect"
 res@vpWidthF      = 0.8   
 res@vpHeightF     = 0.4   
 ;res@vpXF          = 0.05
 ;res@vpYF          = 0.95

  res@gsnMajorLatSpacing = 30              ; change maj lat tm spacing
  res@gsnMajorLonSpacing = 60              ; change maj lon tm spacing
  res@tmXBMinorOn        = False           ; no lon minor tickmarks
  res@tmYLMinorOn        = False
 ; Set resources for tickmarks
  ;res@tmYLLabelsOn               = False            ; turn off lat labels
  ;res@tmXBLabelsOn               = False            ; turn off lon labels

  res@tmYROn                     = False
  ;res@tmYLOn                     = False
  res@tmXTOn                     = False
  ;res@tmXBOn                     = False

; to have a common label bar, all plots should be set to the same interval
; b/c the label bar is drawn by default from the interval of the first plot.
  res@cnLevelSelectionMode =  "ManualLevels"   
  res@cnMinLevelValF       = -1. ; tonnes/ha
  res@cnMaxLevelValF       = 1.  ; 
  res@cnLevelSpacingF      = 0.1 ;       

; change a default color map
  ;---Reading the default colormap into a two dimensional 254 x 4 (RGBA) array
  cmap = read_colormap_file("NEO_div_vegetation_a")
;  cmap = read_colormap_file("ncl_default")
;  ;cmap = cmap(::-1,:) ; reverse color map
;  ;select the middle cmap elements that correspond to 0 in the middle of value range (-0.5~0.5)
;  cmap(110:136,:) = 1 ; use white for zeros; 0 for black, 1 for white
  res@cnMissingValFillColor = "white"  ; missing values (default gray)
  res@cnFillPalette         = cmap     ; use the modified color map; cnFillPalette spans the color map by default

; res@cnLevelSelectionMode = "ExplicitLevels"
; res@cnLevels    = (/-1., -0.75, -0.5, -0.25, -0.001, 0.001, 0.25, 0.5, 0.75, 1./) ; level to divide data into ranges
; ;res@cnFillColors = (/ 0, 54, 2, 3, 4, 5, 28, 20, 14, 0/)
  ;res@cnSpanFillPalette     = True    ; spans cnFillPalette to get color index values for cnFillColors in Explicit mode

; casename = (/ "RCP45~B~(SSP5LU)~N~ - RCP85","RCP85+SAI - RCP85","RCP85+MSB - RCP85","RCP85+CCT - RCP85" /)
; casename = (/ "RCP45 - RCP85","SAI - RCP85","MSB - RCP85","CCT - RCP85" /) 
 casename = (/ "SAI - RCP85","MSB - RCP85","CCT - RCP85", "Emissions reduction" /)
 seq = (/ "(~F22~a~F~) SAI - RCP85", "(~F22~b~F~) MSB - RCP85", "(~F22~c~F~) CCT - RCP85", "(~F22~d~F~) Emissions reduction" /) ; F22 is helvetica bold font

 ;yield difference per case
 ;=========================
  do kk=0,3
    res@gsnLeftString      = seq(kk)
    ;res@gsnRightString     = casename(kk)
    plot1(kk)              = gsn_csm_contour_map_ce(wks1,totyc1(kk+1,:,:),res)
  end do



; create points for box
;************************************************
;set regions: 
 regs = (/"NAM","SEA","EUA","SAM","CAM","SAF"/) 
; latn = (/ 60,  45,  60,     0,  30,   15/)
; lats = (/ 30, -10,  36,   -40,  0,   -30/)
; lonw = (/-130, 80, -10,   -80, -110, -5/)
; lone = (/ -70, 125, 70,   -35, -50,   50/)
;larger box for SEA to include India and SAF
 latn = (/ 60,  55,  60,     0,  30,   15/)
 lats = (/ 30, -10,  36,   -40,  0,   -30/)
 lonw = (/-130, 70, -10,   -80, -110, -12/)
 lone = (/ -70, 130, 70,   -35, -50,   50/)

 nreg = dimsizes(lats);


  resp                  = True
  resp@gsLineDashPattern = 0
  resp@gsLineThicknessF = 1.5
  resp@gsLineColor      = 1
  txres               = True
  txres@txFontHeightF = 0.018             ; Set the font height

  tx = new((/nreg,4/),graphic) ;;---Arrays to hold text annotation
  d1 = new((/nreg,4/),graphic)
  d2 = new((/nreg,4/),graphic)
  d3 = new((/nreg,4/),graphic)
  d4 = new((/nreg,4/),graphic)
  do jj=0,nreg-1
  xpts := (/lonw(jj), lone(jj), lone(jj), lonw(jj), lonw(jj)/)
  ypts := (/latn(jj), latn(jj), lats(jj), lats(jj), latn(jj)/)
  do i = 0,3
  d1(jj,i)=gsn_add_polyline(wks1,plot1(0),xpts(i:i+1),ypts(i:i+1),resp)
  d2(jj,i)=gsn_add_polyline(wks1,plot1(1),xpts(i:i+1),ypts(i:i+1),resp)
  d3(jj,i)=gsn_add_polyline(wks1,plot1(2),xpts(i:i+1),ypts(i:i+1),resp)
  d4(jj,i)=gsn_add_polyline(wks1,plot1(3),xpts(i:i+1),ypts(i:i+1),resp)
  end do
  end do

  ;add label for region boxes
  xlb = lonw + (/-15, 12, -14, -15, -15, 12 /)
  ylb = lats + 8 
  do i = 0,0
   ;tx(:,i) = gsn_add_text(wks1,plot1(i),regs,(lonw+lone)/2,(latn+lats)/2, txres)
   tx(:,i) = gsn_add_text(wks1,plot1(i),regs,xlb,ylb, txres)
  end do

 
 ;===============================
 ;panel plots
  resP                  = True
  ;resP@gsnDraw          = False        ; Turn off draw and frame so
  resP@gsnFrame         = False         ; we can attach some text.
  resP@gsnPanelYWhiteSpacePercent = 5   ;add 5% white space between plots
  ;resP@gsnPanelMainString = "Crop yield change relative to RCP85";
  resP@gsnPanelLabelBar = True ;add a common label to whole panel
  resP@lbLabelFontHeightF = 0.011 ;set label tick font size
  resP@gsnMaximize      = True          ; use full page 
  resP@gsnPaperOrientation = "landscape"; "portrait"
  ;resP@gsnPaperWidthF  = 8.5;11
  ;resP@gsnPaperHeightF = 8.5   
  resP@gsnPaperMargin  = 0             ;minimize margins for SVG plot
  ;resP@gsnPanelBottom   = 0.02         ;add some pace in the bottom (default 0)
  gsn_panel(wks1,plot1,(/2,2/),resP)    ; now draw as on
  
; Draw a text string at the bottom below label bar
  txres               = True
  txres@txFontHeightF = 0.012
  gsn_text_ndc(wks1,"~F8~D~F~Yield (tonnes ha~S~-1~N~)",0.5,0.15,txres) 
  ;function codes: ~F8~D~ is for big delta; ~F~ or ~N~ means switch back to default font
  ;gsn_text_ndc(wks,"~F8~D~",0.5,0.1,txres) ;~F8~D~ is the special code for delta

  frame(wks1)
  ;maximize_output(wks1,resP) 


end
